13 Advanced topic: Porous media wave propagation
We want to consider another multi-physics problem that is motivated from several fields: wave propagation through porous media. Up to now we only considered wave propagation through continuous solids (Chapter 6). But in some cases fluids inside a porous medium make it difficult to describe the medium as purely elastic anymore. One example is the search for hydrocarbons in the subsurface and a second example is the investigation of bones. The theory presented in the following dates back to the seminal work of M. A. Biot in the 1950s.
Consider a representative elementary volume \(\mathit{\Omega}\), which consists of a porous skeleton with volume \(\mathit{\Omega}_\mathrm{s}\) and pores \(\mathit{\Omega}_\mathrm{f}\) that are fully saturated by a fluid (see Figure 13.1) . Note, that this theory can be extended to unsaturated (i.e., pores filled with two or more immiscible fluids) porous media, but this is out of the scope here. The model volume \(\mathit{\Omega}\) has the initial volume \(\mathit{\Omega}_0\) and the current volume \(\mathit{\Omega}_\mathrm{t}\). The Eulerian porosity \(\varphi\) of the macroscopic medium is then defined by
\[ \varphi = \frac{\mathit{\Omega}_\mathrm{f}}{\mathit{\Omega}_\mathrm{t}} = 1 - \frac{\mathit{\Omega}_\mathrm{s}}{\mathit{\Omega}_\mathrm{t}}, \tag{13.1}\]
and the Lagrangian porosity \(\phi\) is defined by
\[ \phi = \frac{\mathit{\Omega}_\mathrm{f}}{\mathit{\Omega}_0} = 1 - \frac{\mathit{\Omega}_\mathrm{s}}{\mathit{\Omega}_0}. \]
To switch between the Lagrangian and Eulerian porosities, the relation
\[ \phi = J \varphi, \quad \text{with} \quad J=\frac{\mathit{\Omega}_\mathrm{t}}{\mathit{\Omega}_0}, \tag{13.2}\]
is used, where \(J\) is the Jacobian of the deformation. Under the approximation of infinitesimal transformation, the Jacobian \(J\) is given by
\[ J \approx 1 + \epsilon, \tag{13.3}\]
where \(\epsilon\) is the linearised volume dilation of the skeleton (i.e. the trace of the linearised strain tensor \(\mathbf{D}\)):
\[ \epsilon = D_{ii} = \mathbf{\nabla} \cdot \mathbf{u}_\mathrm{s}. \tag{13.4}\]
Here, we follow a continuum approach, that means, we do not distinguish at which point in the medium there is a pore and where there is the solid matrix. Each particle point has all the properties and only the porosity defines how much fluid and how much solid there is.
13.1 Balance equations
13.1.1 Mass balance
Recall the mass balance for the solid Equation 3.7 that — for a poroelastic medium — we have to set up for both the fluid and the solid parts separately, i.e.,
\[ \begin{align*} \partial_t \varphi\rho_\mathrm{f} + \nabla \cdot \left( \varphi\rho_\mathrm{f} \mathbf v_\mathrm{f} \right) & = 0, \\ \partial_t \left(1-\varphi\right)\rho_\mathrm{s} + \nabla \cdot \left( \left(1-\varphi\right)\rho_\mathrm{s} \mathbf v_\mathrm{s} \right) & = 0. \end{align*} \]
For the fluid part, we now assume homogeneity on a micro-scale, i.e., \(\nabla \phi\) and \(\nabla\rho_\mathrm{f}\) vanish, and replace the density with the bulk modulus \(K_\mathrm{f}=\rho_\mathrm{f}\frac{\partial p}{\partial \rho_\mathrm{f}}\), where \(p\) is the fluid pressure, to obtain
\[ \partial_t \varphi + \frac{\varphi}{K_\mathrm{f}} \partial_t p + \varphi \nabla \cdot \mathbf v_\mathrm{f} = 0. \tag{13.5}\]
13.1.2 Energy balance
We need to find a relation for the change of porosity \(\partial_t \varphi\). This can be found from the balance of energy. The first law of thermodynamics states that the change in the internal energy \(U\) of a closed system is equal to the amount of heat \(\delta Q\) supplied to the system, plus the amount of work \(\delta W\) done on the system. Mathematically, this law is written as
\[ \mathrm{d}U = \delta Q + \delta W. \]
Setting up the balance of energy for the skeleton, the heat supplied to the system is given by \(\delta Q=T \mathrm{d}\mathcal{S}\), where \(T\) and \(\mathcal{S}\) are the temperature and the entropy of the system, respectively, and the work is composed of the elastic work done by the whole system, \(\delta W_1 = \mathit{\Omega}_0 \mathbf{\sigma} : \mathrm{d}\mathbf{D}\) (where \(\mathbf{\sigma} : \mathrm{d}\mathbf{D} = \sigma_{ij} D_{ij}\)), minus the elastic work done on the fluid, \(\delta W_2 = - p \mathrm{d}\mathit{\Omega}\). In addition, the dissipation, or the change in internal energy \(\mathrm{d}U\) is zero.
This results in
\[ \begin{align*} T\mathrm{d}\mathcal{S} + \mathit{\Omega}_0 \mathbf{\sigma} : \mathrm{d}\mathbf{D} + p \mathrm{d}\mathit{\Omega}_\mathrm{f} & = 0 \\ \Leftrightarrow T\mathrm{d}\mathcal{S}_\mathit{\Omega} + \mathbf{\sigma} : \mathrm{d}\mathbf{D} + p \mathrm{d}\phi & = 0, \end{align*} \]
with \(\mathcal{S}_\mathit{\Omega}\) denoting the specific entropy per unit volume. The term \(T\mathrm{d}\mathcal{S}_\mathit{\Omega}\) may be replaced using the definition of the skeleton free energy \(\mathit{\psi}_\mathit{\Omega}\), which is given by
\[ \mathit{\psi}_\mathit{\Omega} = U_\mathit{\Omega} - T \mathcal{S}_\mathit{\Omega}, \]
where \(U_\mathit{\Omega}\) is the mechanical dissipation of the skeleton per unit volume which is zero for elastic processes. This replacement yields
\[ \mathbf{\sigma} : \mathrm{d}\mathbf{D} + p \mathrm{d}\phi -\mathcal{S}_\mathit{\Omega}\mathrm{d}T - \mathrm{d}\mathit{\psi}_\mathit{\Omega} = 0. \tag{13.6}\]
With the use of the specific Gibbs free energy \(G_\mathrm{s}\), which is given by
\[ G_\mathrm{s} = \psi_\mathrm{s} - p \phi , \]
and assuming that the temperature is constant, we obtain the final energy balance equation for a saturated porous medium:
\[ \mathbf{\sigma} : \mathrm{d}\mathbf{D} - \phi \mathrm{d}p - \mathrm{d}G_\mathrm{s} = 0. \tag{13.7}\]
From equation Equation 13.7, the following state equations can be derived:
\[ G_\mathrm{s} = G_\mathrm{s} \left(D_{ij},p\right); \quad \sigma_{ij}=\frac{\partial G_\mathrm{s}}{\partial D_{ij}}; \quad \phi = -\frac{\partial G_\mathrm{s}}{\partial p}. \tag{13.8}\]
For a more detailed analysis of these state equations, see (Coussy 2003).
Calculating the total derivatives of the state variables and using the state equations Equation 13.8 recovers the following expressions:
\[ \begin{align*} \mathrm{d}\sigma_{ij} & = \frac{\partial \sigma_{ij}}{\partial D_{kl}} \mathrm{d}D_{kl} + \frac{\partial \sigma_{ij}}{\partial p} \mathrm{d}p \\* \Leftrightarrow \mathrm{d}\sigma_{ij} & = \frac{\partial^2 G_\mathrm{s}}{\partial D_{ij}\partial D_{kl}} \mathrm{d}D_{kl} + \frac{\partial^2 G_\mathrm{s}}{\partial p \partial D_{ij}} \mathrm{d}p, \\ \mathrm{d}p & = \frac{\partial \phi}{\partial D_{ij}} \mathrm{d}D_{ij} + \frac{\partial \phi}{\partial p} \mathrm{d}p \\* \Leftrightarrow \mathrm{d}\phi & = - \frac{\partial^2 G_\mathrm{s}}{\partial D_{ij}\partial p} \mathrm{d}D_{ij} - \frac{\partial^2 G_\mathrm{s}}{\partial p^2} \mathrm{d}p. \end{align*} \]
This leads to the incremental state equations of thermoporoelasticity:
\[ \begin{align*} \mathrm{d}\sigma_{ij} & = E_{ijkl} \mathrm{d}D_{kl} - b_{ij} \mathrm{d}p \tag*{🟦}, \\ \mathrm{d}\phi & = b_{ij} \mathrm{d}D_{ij} + \frac{1}{N} \mathrm{d}p \tag*{🟧}, \end{align*} \tag{13.9}\]
where \(E_{ijkl} = \frac{\partial^2 G_\mathrm{s}}{\partial \varepsilon_{ij} \partial \varepsilon_{kl}}\) is the (drained) elastic tensor, \(b_{ij}=-\frac{\partial^2 G_\mathrm{s}}{\partial \varepsilon_{ij} \partial p}\) is the Biot effective stress coefficient tensor , and \(\frac{1}{N} = -\frac{\partial^2 G_\mathrm{s}}{\partial p^2}\) is the inverse of Biot’s modulus (Coussy 2003). These properties are often also called tangent properties because they define the derivatives of the state equation, ergo the slope of the tangent, and are set constant as usually done in linear thermoporoelasticity. While most properties are familiar or sufficiently described by their name, it is worth having a closer look at the inverse of Biot’s modulus. This property links the pressure variation \(\mathrm{d}p\) and the porosity variation \(\mathrm{d}\phi\) if strain is held constant (\(\mathrm{d}\mathbf{D}=\mathbf 0\)) via the expression \(\mathrm{d}\phi=N^{-1} \mathrm{d}p\). It, therefore, tells how easy it is to enlarge the pores when the pressure is increased.
We can integrate these equations and obtain: \[ \begin{align*} \sigma_{ij} - \sigma^0_{ij} & = E_{ijkl} \left(D_{kl} - D^0_{kl}\right) - b_{ij} \left(p-p^0\right) \\ \phi-\phi^0 & = b_{ij} \left(D_{kl} - D^0_{kl}\right) + \frac{1}{N} \left(p-p^0\right). \end{align*} \]
Afterwards we can take the time derivative and get
\[ \begin{align*} \partial_t \sigma_{ij} & = E_{ijkl} \partial_t D_{kl} - b_{ij} \partial_t p \tag*{🟦}, \\ \partial_t \phi & = b_{ij} \partial_t D_{ij} + \frac{1}{N} \partial_t p \tag*{🟧}, \end{align*} \tag{13.10}\]
Equation Equation 13.10 \(🟧\) can also be written for the Eulerian porosity. With \(\phi = J \varphi\) (Equation 13.2) and \(J \approx 1 + \epsilon\) (Equation 13.3),
\[ \partial_t \phi = \partial_t \left( J \varphi \right) = J \partial_t \varphi + \varphi \partial_t \epsilon \simeq \partial_t \varphi + \varphi \partial_t \epsilon , \]
where the last step is valid for small deformations, where \(J \approx 1\). Hence, equation Equation 13.10 \(🟧\) using equation Equation 13.4 becomes
\[ \partial_t \varphi = -\varphi \mathbf{\nabla} \cdot \mathbf{v}_\mathrm{s} + b_{ij} \partial_t D_{ij} + \frac{1}{N} \partial_t p. \tag{13.11}\]
This equation is named incremental state equation for porosity and provides a relation for the temporal change of the porosity.
13.1.3 Momentum balance
Recall the momentum balance Equation 3.9,
\[ \rho \frac{D}{Dt} \mathbf v = \nabla \cdot \mathbf \sigma + \rho \mathbf b. \tag{13.12}\]
This can be used to describe a porous medium on the macro scale, too. However, we also have to consider the momentum balance on the micro scale, i.e., the interactions between the fluid and the solid matrix. Hence, within the different volumes of the solid and fluid, the momentum balance can be written as
\[ \begin{align*} \int\limits_{\mathit{\Omega}_\mathrm{s}} \rho_\mathrm{s} \frac{D}{Dt} \mathbf v_\mathrm{s} \mathrm{d}\mathit{\Omega}_\mathrm{s} & = \int\limits_{\mathit{\Omega}\mathrm{s}} \nabla \cdot \mathbf \sigma_\mathrm{s} \mathrm{d}\mathit{\Omega}_\mathrm{s} + \int\limits_{\mathit{\Omega}\mathrm{s}} \rho_\mathrm{s} \mathbf b \mathrm{d}\mathit{\Omega}_\mathrm{s}, \\ \int\limits_{\mathit{\Omega}_\mathrm{f}} \rho_\mathrm{f} \frac{D}{Dt} \mathbf v_\mathrm{f} \mathrm{d}\mathit{\Omega}_\mathrm{f} & = \int\limits_{\mathit{\Omega}_\mathrm{f}} \nabla \cdot \mathbf \sigma_\mathrm{f} \mathrm{d}\mathit{\Omega}_\mathrm{f} + \int\limits_{\mathit{\Omega}_\mathrm{f}} \rho_\mathrm{f} \mathbf b \mathrm{d}\mathit{\Omega}_\mathrm{f}. \end{align*} \]
We can add the two equations to consider the interactions between the fluid and the solid. Furthermore, we consider that the quantities in the whole medium \(\mathit{\Omega}\) can be expressed as the quantities in the partial volumes \(\mathit{\Omega}_\mathrm{s}\) and \(\mathit{\Omega}_\mathrm{f}\) times the corresponding volume fraction, i.e. we average the over the volume \(\mathit{\Omega}\):
\[ \begin{gathered} \int\limits_{\mathit{\Omega}} \left(1-\varphi\right) \rho_\mathrm{s} \mathbf a_\mathrm{s} \mathrm{d}\mathit{\Omega} + \int\limits_{\mathit{\Omega}} \varphi \rho_\mathrm{f} \mathbf a_\mathrm{f} \mathrm{d}\mathit{\Omega} = \\ \int\limits_{\mathit{\Omega}} \left(1-\varphi\right) \rho_\mathrm{s} \mathbf{b} \mathrm{d}\mathit{\Omega} + \int\limits_{\mathit{\Omega}} \varphi \rho_\mathrm{f} \mathbf{b} \mathrm{d}\mathit{\Omega} + \int\limits_{\mathit{\Omega}} \mathbf{\nabla} \cdot \left[ \left(1-\varphi\right) \mathbf{\sigma}_\mathrm{s} \right] \mathrm{d}\mathit{\Omega} + \int\limits_{\mathit{\Omega}} \mathbf{\nabla} \cdot \left( \varphi \mathbf{\sigma}_\mathrm{f} \right) \mathrm{d}\mathit{\Omega}, \end{gathered} \tag{13.13}\]
where \(\mathrm{a} = \frac{D}{Dt} \mathbf v\) is the acceleration. Note that this acceleration may be recognised as acceleration (in the Eulerian view) or particle-acceleration (in the Lagrangian view). It represents a physical acceleration. Therefore, the particle derivative also serves as a link between the Eulerian and Lagrangian descriptions. Because we can choose \(\mathit{\Omega}\) arbitrarily, we can write
\[ \left(1-\varphi\right) \rho_\mathrm{s} \mathbf a_\mathrm{s} + \varphi \rho_\mathrm{f} \mathbf a_\mathrm{f} = \left[\left(1-\varphi\right) \rho_\mathrm{s} + \varphi \rho_\mathrm{f} \right] \mathbf{b} + \mathbf{\nabla} \cdot \left[ \left(1-\varphi\right) \mathbf{\sigma}_\mathrm{s} + \varphi \mathbf{\sigma}_\mathrm{f} \right]. \tag{13.14}\]
Upon comparison with the microscopic momentum balance and the macroscopic momentum balance Equation 13.12, we can infer the stress partition theorem
\[ \mathbf{\sigma} = \left( 1-\varphi \right) \mathbf{\sigma}_\mathrm{s} + \varphi \mathbf{\sigma}_\mathrm{f} = \left( 1-\varphi \right) \mathbf{\sigma}_\mathrm{s} - \varphi p \tag{13.15}\]
and the definition of the bulk density
\[ \rho = \left( 1 - \varphi \right) \rho_\mathrm{s} + \varphi \rho_\mathrm{f}. \tag{13.16}\]
We, hence, can write the momentum balance for a poroelastic medium as
\[ \mathbf{\nabla} \cdot \mathbf{\sigma} - \left( 1-\varphi \right) \rho_\mathrm{s} \partial_t \mathbf{v}_\mathrm{s} - \varphi \rho_\mathrm{f} \partial_t \mathbf{v}_\mathrm{f} + \rho \mathbf{b} = \mathbf{0} . \tag{13.17}\]
The body force density \(\rho \mathbf b\) can usually be reduced to a source term \(\rho \mathbf s\), because other body forces, which in particular include gravity, can usually be neglected for wave propagation problems. This resulting equation is often called Biot’s dynamic equation.
13.2 Constitutive equations
13.2.1 Generalised Hooke’s law
The generalised Hooke’s law can be derived from the balance of mass and the balance of energy. Inserting the incremental state equation for porosity Equation 13.11 into the balance of mass for the fluid Equation 13.5 yields
\[ b_{ij} \partial_t D_{ij} + \frac{1}{N} \partial_t p + \frac{\varphi}{K_\mathrm{f}} \partial_t p + \varphi \nabla \cdot \mathbf v_\mathrm{f} - \varphi \mathbf{\nabla} \cdot \mathbf{v}_\mathrm{s} = 0 , \tag{13.18}\]
which we can also write together with the time-derivative of the first incremental state equation of poroelasticity Equation 13.10 \(🟦\) as
\[ \begin{align*} \partial_t \sigma_{ij} & = E_{ijkl} \partial_t D_{kl} - b_{ij} \partial_t p, \tag*{🟦} \\ \partial_t p & = M \left[ \varphi \left( \mathbf{\nabla} \cdot \mathbf{v}_\mathrm{s} - \nabla \cdot \mathbf v_\mathrm{f} \right) - b_{ij} \partial_t D_{ij} \right] \tag*{🟧} , \end{align*} \tag{13.19}\]
with
\[ \frac{1}{M} = \frac{1}{N} + \frac{\varphi}{K_\mathrm{f}} . \tag{13.20}\]
Inserting equation Equation 13.19 \(🟧\) into equation Equation 13.19 \(🟦\) gives an expression that uses the kinematic quantities, which are the particle velocities \(\mathbf{v}_\mathrm{s}\) and \(\mathbf{v}_\mathrm{f}\), instead of the dynamic quantity \(p\):
\[ \partial_t \sigma_{ij} = E^\mathrm{u}_{ijkl} \partial_t D_{kl} - \varphi M b_{ij} \left( \mathbf{\nabla} \cdot \mathbf{v}_\mathrm{s} - \nabla \cdot \mathbf v_\mathrm{f} \right) \tag{13.21}\]
where
\[ E^\mathrm{u}_{ijkl} = E_{ijkl} + M b_{ij}b_{kl} \tag{13.22}\]
is the undrained elastic tensor.
13.2.2 Fluid transport
Finally, we need an equation that describes the fluid flow inside the pores. Recall Darcy’s law as derived previously Equation 8.17: q \[ \nabla p = \rho \mathbf g - \varphi \frac{\mu}{\kappa} \mathbf v_\mathrm{f}. \]
However, in a poroelastic medium we have to write this in terms of the relative velocity of the fluid with respect to the solid frame. This extended Darcy’s law that describes viscosity-dominated Hagen-Poiseuille flow reads
\[ \nabla p = - \rho_\mathrm{f} \left[ T \partial_t \mathbf v_\mathrm{f} + \left(1-T\right) \partial_t \mathbf v_\mathrm{s} \right] - \varphi \frac{\mu}{\kappa} \left( \mathbf v_\mathrm{f} - \mathbf v_\mathrm{s} \right), \tag{13.23}\]
where \(T\) is the tortuosity, defined as the ratio of the effective length of a flow path from a point A to another point B and the direct distance between these points. It is usually determined experimentally or the relation obtained by (Berryman 1980),
\[ T = 1 - \mathfrak{r} \left( 1 - \frac{1}{\varphi} \right) , \tag{13.24}\]
is used, where \(\mathfrak{r}\) is a factor that depends on the shape of the pores. In general, \(\mathfrak{r}\) ranges between \(0\) and \(1\) for arbitrary ellipsoids and, for spheres, \(\mathfrak{r}=\frac{1}{2}\). The tortuosity is equal to \(1\) for cylindrical pores with axes parallel to the macroscopic pressure gradient and approaches infinity for a very low porosity.
13.3 Wave equation
Finally, the governing equations for wave propagation in saturated porous media are the generalised Hooke’s law given by equations Equation 13.21 and Equation 13.19 \(🟧\), the extended Darcy’s law as previously presented in equation Equation 13.23 and Biot’s dynamic equation as given in equation Equation 13.17. They are, using the index notation, written as
\[ \begin{align*} \partial_t \sigma_{ij} = & E^\mathrm{u}_{ijkl} \partial_t D_{kl} - \varphi M b_{ij} \left( \partial_k v_{\mathrm{s}k} - \partial_k v_{\mathrm{f}k}\right), \tag*{🟦} \\ % \partial_t p = & M \left[ \varphi \left( \partial_i v_{\mathrm{s}i} - \partial_i v_{\mathrm{f}i} \right) - b_{ij} \partial_t D_{ij} \right], \tag*{🟧} \\ % \partial_i p = & - \rho_\mathrm{f} \left[ T \partial_t v_{\mathrm{f}i} + \left(1-T\right) \partial_t v_{\mathrm{s}i} \right] - \varphi \frac{\mu}{\kappa} \left( v_{\mathrm{f}i} - v_{\mathrm{s}i} \right), \tag*{🟨} \\ % \partial_j \sigma_{ij} = & \left( 1-\varphi \right) \rho_\mathrm{s} \partial_t v_{\mathrm{s}i} + \varphi \rho_\mathrm{f} \partial_t v_{\mathrm{f}i} - \rho s_i \tag*{🟩} . \end{align*} \tag{13.25}\]
We can rearrange the equations and combine Equation 13.25 \(🟦\) and \(🟧\) to get the final system of partial differential equations in a velocity-stress formulation:
\[ \begin{align*} \partial_t \sigma_{ij} & - E^\mathrm{u}_{ijkl} \partial_k v_{\mathrm{s}l} + \varphi M b_{ij} \partial_k v_{\mathrm{s}k} - \varphi M b_{ij} \partial_k v_{\mathrm{f}k} = 0, \tag*{🟦} \\ % \partial_t p & + M b_{ij} \partial_i v_{\mathrm{s}j} - \varphi M \partial_i v_{\mathrm{s}i} + \varphi M \partial_i v_{\mathrm{f}i} = 0, \tag*{🟧} \\ % \partial_t v_{\mathrm{s}i} & - \frac{1}{\rho^*} \partial_j \sigma_{ij} - \frac{\varphi}{T\rho^*} \partial_i p = \frac{\varphi^2 \mu}{T \rho^* \kappa} \left( v_{\mathrm{f}i} - v_{\mathrm{s}i} \right) + \frac{\rho}{\rho^*} s_i, \tag*{🟨} \\ % \partial_t v_{\mathrm{f}i} & + \frac{1-T}{T\rho^*} \partial_j \sigma_{ij} + \frac{\left(1-\varphi\right)\rho_\mathrm{s}}{T\rho_\mathrm{f}\rho^*} \partial_i p = - \frac{\left(1-\varphi\right)\rho_\mathrm{s}}{T\rho_\mathrm{f}\rho^*}\frac{\varphi\mu}{\kappa} \left( v_{\mathrm{f}i} - v_{\mathrm{s}i} \right) - \frac{\left(1-T\right)\rho}{T\rho^*} s_i \tag*{🟩} , \end{align*} \tag{13.26}\]
with the effective density
\[ \rho^* = \left(1-\varphi\right) \rho_\mathrm{s} + \left(1-\frac{1}{T} \right) \varphi \rho_\mathrm{f}. \]
If the porous medium is isotropic, the Biot tensor is given by
\[ b_{ij} = b \delta_{ij} = \left( 1-\frac{K}{K_\mathrm{s}}\right) \delta_{ij}, \]
where \(K\) is the skeleton bulk modulus and \(K_\mathrm{s}\) is the bulk modulus of the material that builds up the skeleton. In addition, Biot’s modulus can be written as
\[ \frac{1}{N} = \frac{b-\varphi}{K_\mathrm{s}}. \]
The isotropic elastic modulus \(E_{ijkl}\) only consists of two independent parameters (e.g., the Lamé parameters \(\lambda\) and \(\mu\)) according to
\[ E_{ijkl} = E_{IJ} = \left( \begin{matrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{matrix} \right) , \]
where the Voigt notation is used.
In contrast to the elastic wave equation in velocity-stress notation
\[ \begin{align*} \partial_t \sigma_{ij} - E_{ijkl} \partial_k v_l & = 0 \tag*{🟦} \\ \partial_t v_i - \frac{1}{\rho} \partial_j \sigma_{ij} & = 0 \tag*{🟧}, \end{align*} \]{eq-waveprop-hooke_new}
we now have four additional equations describing the fluid velocity and the fluid pressure. Additionally, we have terms with the velocity fields themselves (i.e., no time or space derivative). These terms make the equation an advection-reaction equation (in constrast to the simple advection equation for elastic wave propagation) and are hence called reaction terms. These reaction terms cause damping of the waves.
13.3.1 Special case of elastic wave propagation
If this theory that describes wave propagation in porous media is well derived, it has to be possible to deduce theories that describe limiting cases of this theory. Hence, it has to be possible to deduce the wave equation for elastic media by calculating the limit of vanishing porosity. First, the effective density, becomes
\[ \lim\limits_{\varphi \rightarrow 0} \rho^* = \lim\limits_{\varphi \rightarrow 0} \left[ \left(1-\varphi\right) \rho_\mathrm{s} + \left(1-\frac{1}{T} \right) \varphi \rho_\mathrm{f} \right] = \rho_\mathrm{s}, \]
and the bulk density becomes
\[ \lim\limits_{\varphi \rightarrow 0} \rho = \lim\limits_{\varphi \rightarrow 0} \left[ \left( 1 - \varphi \right) \rho_\mathrm{s} + \varphi \rho_\mathrm{f} \right] = \rho_\mathrm{s}. \]
The tortuosity implicitly depends on the porosity as well. With decreasing porosity, the flow path should become more and more tortuous, leading to an infinite tortuosity as the porosity approaches zero. Using equation Equation 13.24, this can be demonstrated:
\[ \lim\limits_{\varphi \rightarrow 0} T = \lim\limits_{\varphi \rightarrow 0} \left[ 1 - r \left( 1 - \frac{1}{\varphi} \right) \right] = \infty . \]
Hence, \(\lim\limits_{\varphi \rightarrow 0} T^{-1}=0\). Next, it is necessary to have a closer look at Biot’s effective stress coefficient tensor \(b_{ij}\). This parameter also intrinsically depends on the porosity. The skeleton bulk modulus \(K\) and the matrix bulk modulus \(K_\mathrm{s}\) must become the same as the porosity approaches zero. Without any pores, the skeleton is the same as the material that builds the skeleton. Hence, also the bulk properties should be the same. Therefore,
\[ b = 1- \frac{K}{K_\mathrm{s}} = 1- \frac{K_\mathrm{s}}{K_\mathrm{s}} = 0 . \]
Furthermore, it follows that \(E^\mathrm{u}_{ijkl} = E_{ijkl}\).
With that, the system of partial differential equation reduces to
\[ \begin{align*} \partial_t \sigma_{ij} & - E_{ijkl} \partial_k v_{\mathrm{s}l} = 0, \tag*{🟦} \\ % \partial_t p & = 0, \tag*{🟧} \\ % \partial_t v_{\mathrm{s}i} & - \frac{1}{\rho_\mathrm{s}} \partial_j \sigma_{ij} = s_i \tag*{🟨}, \\ % \partial_t v_{\mathrm{f}i} & - \frac{1}{\rho_\mathrm{s}} \partial_j \sigma_{ij} = s_i \tag*{🟩} . \end{align*} \tag{13.27}\]
The reaction terms become zero, which makes sense since it accounts for energy dissipation and viscodynamic effects. These are not present in elastic wave equations. These equations coincide with the common elastic wave equations while the velocity fields are just becoming the same and describe the same field and the pressure field vanished.
13.4 Propagation of Harmonic Plane Waves
The system of partial differential equations Equation 13.26 (\(🟦, 🟧, 🟨, 🟩\)) can be written in the matrix-vector form
\[ \partial_t Q_p + A_{pq} \partial_x Q_q + B_{pq} \partial_y Q_q + C_{pq} \partial_z Q_q = E_{pq} Q_q + R_{pr} s_r, \tag{13.28}\]
where \(\mathbf{Q}=\left(\sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{xy}, \sigma_{yz}, \sigma_{xz}, v_{\mathrm{s} x}, v_{\mathrm{s} y}, v_{\mathrm{s} z}, p, v_{\mathrm{f} x}, v_{\mathrm{f} y}, v_{\mathrm{f} z} \right)^\mathrm{T}\) is the solution vector containing the \(13\) unknown variables. \(\mathbf{A}\), \(\mathbf{B}\) and \(\mathbf{C}\) are the space-dependent matrices of size \(13\times13\) that contain the material properties. \(\mathbf{E}\) is the reaction term, which accounts for energy dissipation and viscodynamic effects.
The propagation of harmonic plane waves is considered by choosing the following ansatz:
\[ Q_\mathrm{q} \left(\mathbf{x},t\right) = a_\mathrm{q} \exp\left[ \mathrm{i} \left( \mathbf{k} \cdot \mathbf{x} - \omega t + \phi_0 \right) \right] , \]
where \(\mathbf{k}\) is the complex wavenumber vector, \(\omega\) is the angular frequency, \(\phi_0\) is the initial phase shift and \(a_q\) is the amplitude. Substituting this ansatz into the previously derived wave equation (Equation 13.28) yields
\[ \left( -\omega \delta_{pq} + k_x A_{pq} + k_y B_{pq} + k_z C_{pq} + \mathrm{i} E_{pq} \right) a_q \exp\left[ \mathrm{i} \left( \mathbf{k} \cdot \mathbf{x} - \omega t + \phi_0 \right) \right] = 0 . \]
Since \(\exp\left[ \mathrm{i} \left( \mathbf{k} \cdot \mathbf{x} - \omega t + \phi_0 \right) \right] \neq 0\), the equation reduces to
\[ \left( -\omega \delta_{pq} + k_x A_{pq} + k_y B_{pq} + k_z C_{pq} + \mathrm{i} E_{pq} \right) a_q = 0 . \]
For simplicity, it is assumed that the plane waves propagate in \(x\)-direction only. Therefore, the equation reduces further to
\[ \left( -\omega \delta_{pq} + k_x A_{pq} + \mathrm{i} E_{pq} \right) a_q = 0 . \tag{13.29}\]
This is basically a simple eigenvalue problem. To obtain the dispersion relation for the harmonic waves, this equation has to be solved for the wavenumber \(k_x\) as a function of frequency \(\omega\):
\[ \det \left( -\omega \delta_{pq} + k_x A_{pq} + \mathrm{i} E_{pq} \right) = 0 . \tag{13.30}\]
The resulting polynomial in \(k_x\) has a degree of \(13\) but, because the rank of \(\mathbf{A}\) is eight, five solutions are equal to zero. The remaining polynomial with a degree of eight is even and, therefore, has only four different values with the same norm. The sign just represents the direction of propagation. The remaining four values are the four wavenumbers that correspond to four waves: two P-waves and two S-waves. The two P-waves are the ‘normal’ or fast P-wave (P1-wave) and an additional slow P2-wave.x The two S-waves are the ‘normal’ SH- and SV-waves, i.e., horizontally and vertically polarized shear waves, whose displacement vectors are orthogonal to each other.
13.4.1 Dispersion Relations
If waves propagate at different velocities for varying frequencies or wavelengths, this behaviour is called dispersion. This causes, for example, that a wave packet of a range of frequencies can spread out in space while propagating. The velocity \(c\) of a wave is, obviously, a function of the frequency \(f\):
\[ c\left(f\right) = \mathfrak{F}\left(f\right) , \]
where \(\mathfrak{F}\left(f\right)\) is the dispersion relation. Of course, this relation can also be expressed in terms of the angular frequency \(\omega\). To calculate dispersion relations, Equation 13.30 is solved and the wave velocities are calculated as
\[ c\left(\omega\right) = \frac{\Re\left(\omega\right)}{\Re\left(k_x\left(\omega\right)\right)} . \tag{13.31}\]
13.4.2 Attenuation
Attenuation is a term that is used to describe energy dissipation by all kinds of physical mechanisms that cause a wave amplitude reduction within a material. This does not include geometric wavefront spreading. There are several kinds of attenuation parameters available, but most often attenuation is conventionally described by quality factors, \(Q\). These quality factors are attributes of the materials and defined analogously to their counterparts in electrical circuit theory as
\[ \frac{2\pi}{Q} = \frac{\Delta E}{E} , \]
where \(E\) is the energy stored in an oscillating system and \(\Delta E\) is the dissipated energy per cycle. In practice, attenuation is often presented by the inverse quality factor \(Q^{-1}\) or by using a logarithmic scale. It is calculated from the solutions of equation Equation 13.30 by
\[ Q^{-1} = 2\left|\frac{\Im\left(k_x\left(\omega\right)\right)}{\Re\left(k_x\left(\omega\right)\right)}\right| . \tag{13.32}\]
Small inverse quality factors \(Q^{-1}\) mean always low attenuation while high inverse quality factors mean high attenuation. This is independent of the type of the considered quality factor. For \(Q^{-1}<2\), a system is said to be underdamped. The waves are still oscillating but the amplitude is reducing. For \(Q^{-1}>2\), a system is said to be overdamped. The system is not oscillating at all but exponentially approaches a steady state. Hence, there is no wave travelling any more but it is a diffusion process that occurs. A critical damping is reached for \(Q^{-1}=2\). In this case, there is no oscillation as in the overdamped system but, in contrast, it approaches the steady state in the quickest possible way.
13.4.3 Phase Shift and Amplitudes
In addition to the wavenumber \(k_x\), the (eigen-)vectors \(\mathbf{a}\) with the amplitudes \(a_q\) in equation Equation 13.29 can be investigated as well. They reveal information about the phase shift between the different components of the material, that is between the fluids and the solid part, and about the amplitude of the waves in the different components. The vectors \(\mathbf{a}\) have the same structure as the solution vector \(\mathbf{Q}\) where each component is a complex number, \(a_q = \hat{a}_q \exp{\left(\mathrm{i}\phi_q\right)}\), with the amplitude \(\hat{a}_q\) and the phase shift \(\phi_q\). It is sufficient to consider only the components of \(\mathbf{a}\) that contain the particle velocities. Since these vectors are only unique except a scalar factor, the amplitudes are normalized using \(\sqrt{\hat{a}_\mathrm{s}^2 + \hat{a}_\mathrm{f}^2}=1\) and the phase angle is conveniently calculated relative to the solid as \(\phi_{\mathrm{rel.}} = \phi_\mathrm{f} - \phi_\mathrm{s}\).
13.4.4 Example
As an example the dispersion curves for a Bentheim Sandstone filled with water are presented in Figure 13.2, Figure 13.3 and Figure 13.4. The physical properties used to calculate these curves are provided in Table 13.1.
| Parameter | Symbol | Unit | Value |
|---|---|---|---|
| Density of the solid | \(\rho_\mathrm{s}\) | kg m−3 | 2622.65 |
| Drained bulk modulus of the skeleton | \(K\) | GPa | 7.63 |
| Shear modulus | \(\mu\) | GPa | 7.18 |
| Bulk modulus of the solid grains | \(K_\mathrm{s}\) | GPa | 36.4 |
| Porosity | \(\varphi\) | - | 0.2372 |
| Permeability | \(k\) | m2 | 1.5 × 10−13 |
| Tortuosity | \(T\) | - | 2.608 |
| Density of the fluid | \(\rho_\mathrm{f}\) | kg m−3 | 998.2 |
| Bulk modulus of the fluid | \(K_\mathrm{f}\) | MPa | 2.08 × 109 |
| Viscosity of the fluid | \(\eta_\mathrm{f}\) | Pa s | 1.00 × 10−3 |
Let’s first consider the “normal” P-wave (Figure 13.2) and the S-wave (Figure 13.3). It can be seen that the waves indeed show a frequency dependence, which means that they are dispersive. In general, the propagation velocity is higher at higher frequencies and lower at lower frequencies. Only in a quite narrow transition zone, the waves are damped (\(Q^{-1} \neq 0\)). From the phase shift, we can see the reason of this damping. In general, the solid and fluid phases are moving in phase, but in the transition frequency range, there is out-of-phase movement that causes friction and, hence, energy losses. It is also remarkable that the amplitude of the waves is the same in the low-frequency range and that the waves in the solid part are dominant in the high-frequency range. The additional slow P2-wave (Figure 13.4) is characterized by out-of-phase movement. In general, the solid and liquid parts are moving exactly in opposite direction. This is the reason, why the wave is critically damped (\(Q^{-1} = 2\)) at low frequencies. The amplitude is highly dominating in the fluid. At low frequencies the P2-wave is a diffusive process within the fluid part and not a real wave, which is also the reason for the low velocity. At higher frequencies the P2-wave becomes a normal wave.